clear;
close all;
%% 定义信道参数
global   B c
M_mod = 4;                         % 16QAM
k = log2(M_mod);                    % Bits/symbol
M = 512;
Mzp = 180;
Mp = 0;
B = 5000;                      % 带宽 Hz
scs = 10;                      % 子载波间隔 Hz
ofdmSym = 16;                   % ofdm信号数量 /帧 
codeRate = 1/2;                 % LDPC码率
maxIterations = 20;             % 最大迭代次数LDPC decoder
totalBits = 1e6;                % The approx. total no. of bits simulated
repeats = 1;                    % 重复次数
fc = 12500;
c = 1500;                                                     % speed of light (m/s)
eng_sqrt = (M_mod==2)+(M_mod~=2)*sqrt((M_mod-1)/6*(2^2));
SNR_dB = 5:5:25;
SNR = 10.^(SNR_dB/10);
noise_var_sqrt = sqrt(1./SNR);
sigma_2 = abs(sqrt(eng_sqrt)*noise_var_sqrt).^2;                   % 噪声功率
numSC = M; % 计算载波数量
numDC = (numSC - Mzp-Mp);             % 计算数据长度
% 初始化高斯信道和误码率存储（设置为0）对象
awgnChannel = comm.AWGNChannel('NoiseMethod','Variance', 'VarianceSource','Input port');
errorRate = comm.ErrorRate('ResetInputPort',true);
errorRate1 = comm.ErrorRate('ResetInputPort',true);
errorRate2 = comm.ErrorRate('ResetInputPort',true);
% 初始化LDPC对象
parityCheck_matrix = dvbs2ldpc(codeRate);               % 创建 DVB-S.2 校验矩阵
ldpcEncoder = comm.LDPCEncoder(parityCheck_matrix);     % 编码对象
ldpcDecoder = comm.LDPCDecoder(parityCheck_matrix);     % 解码对象
ldpcDecoder.MaximumIterationCount = maxIterations;      % 分配最大迭代次数
noCodedbits = size(parityCheck_matrix,2);               % 码字比特长度
% 创建误码率储存对象
berOFDM = zeros(length(SNR),3); 
berCOFDM = zeros(length(SNR),3);
berOTFS = zeros(length(SNR),3);
berCOTFS_MRC = zeros(length(SNR),3);
berCOTFS_LMMSE = zeros(length(SNR),3);
berchannal = zeros(length(SNR),3);
errorStats_coded1 = zeros(1,3);
errorStats_coded2 = zeros(1,3);
errorStats_coded3 = zeros(1,3);
errorStats_coded4 = zeros(1,3);
errorStats_coded5 = zeros(1,3);
errorStats_coded6 = zeros(1,3);
N = ofdmSym;
Fn = dftmtx(N);       % DFT matrix
Fn = Fn./norm(Fn);    
Md = M-Mzp-Mp;
data_grid = zeros(M,N);
data_grid(1:Md,:) = 1;
delta_f = scs;
T = 1/delta_f;

max_speed=20;  % km/hr
[chan_coef,delay_taps,Doppler_taps,taps]=Gen_delay_Doppler_channel_parameters(N,M,fc,delta_f,T,max_speed,1);
L_set=unique(delay_taps);
[G,gs]=Gen_time_domain_channel(N,M,taps,delay_taps,Doppler_taps,chan_coef);
[~,K_ml]=Gen_DD_channel_kml(N,M,L_set,gs); 
[H_tf]=Gen_time_frequency_channel(N,M,gs,L_set);

for repetition=1:repeats                
    %生成数据
    [dataIn, dataBits_in, codedData_in, packetSize, numPackets, numCB] = dataGen(k,numDC,ofdmSym,totalBits,codeRate,ldpcEncoder);
     data_na = reshape(dataIn',k*numDC*ofdmSym,packetSize).';
%%                       OTFS
    % 不同EbNo循环
    for m = 1:length(SNR)
        % 不同包循环
        for j = 1:numPackets
            rxframeBuffer1 = [];                 % Initialise matrix 
            rxframeBuffer2 = [];  
            % 不同帧循环
            for u = 1:packetSize
                %一次发送一个帧\               
             data=qammod(reshape(data_na(u,:)',k,Md*N), M_mod,'gray','InputType','bit');
             X = Gen_2D_data_grid(N,M,data,data_grid);     
        % OTFS 
        X_tilda=X*Fn';
        s = reshape(X_tilda,N*M,1);
         % channel output
        r=zeros(N*M,1);
        noise= sqrt(sigma_2(m)/2)*(randn(size(s)) + 1i*randn(size(s)));
        l_max=max(delay_taps);
        for q=1:N*M
            for l=(L_set+1)
                if(q>=l)
                    r(q)=r(q)+gs(l,q)*s(q-l+1);
                end
            end
        end
        r=r+noise;
         % OTFS demodulation
        Y_tilda=reshape(r,M,N);
        Y = Y_tilda*Fn;
        y=reshape(Y.',N*M,1);
            % detector          
            n_ite_MRC = 20;               
            omega = 1;           
            decision = 1;      
            init_estimate = 1;  
            [est_info_bits_LMMSE,~] = LMMSE_detector(N,M,M_mod,sigma_2(m),data_grid,r,gs,L_set); 
            [est_info_bits_MRC_es,~,~]= MRC_detector(N,M,Md,M_mod,sigma_2(m),data_grid,r,H_tf,L_set,omega,decision,init_estimate,n_ite_MRC,K_ml,Y);
            rxframeBuffer1 = [rxframeBuffer1;est_info_bits_MRC_es];     % 储存帧信息
            rxframeBuffer2 = [rxframeBuffer2;est_info_bits_LMMSE];  
            end
            dataOut1 = rxframeBuffer1;
            dataOut2 = rxframeBuffer2;           
            codedData_out1 = randdeintrlv(dataOut1,4831);                      %去交织
            codedData_out1(numel(codedData_in)+1:end) = [];                   % 删除填充比特            
            codedData_out2 = randdeintrlv(dataOut2,4831);                      %去交织
            codedData_out2(numel(codedData_in)+1:end) = [];                   % 删除填充比特          
            dataBits_out1 = [];  
            dataBits_out2 = [];           
            dataOut_buffer1 = codedData_out1;
            dataOut_buffer2 = codedData_out2;
            for q = 1:numCB
                dataBits_out1 = [dataBits_out1;ldpcDecoder(1-2*dataOut_buffer1(1:noCodedbits))]; %解码
                dataOut_buffer1(1:noCodedbits) = [];                     
                dataBits_out2 = [dataBits_out2;ldpcDecoder(1-2*dataOut_buffer2(1:noCodedbits))]; %解码
                dataOut_buffer2(1:noCodedbits) = [];  
            end
            errorStats_coded1 = errorRate1(dataBits_in,dataBits_out1,0);     % 计算错误信息
            errorStats_coded2 = errorRate2(dataBits_in,dataBits_out2,0);            
        end
        berCOTFS_MRC(m,:) =berCOTFS_MRC(m,:)+ errorStats_coded1.';                                  % 保存未编码错误信息
        berCOTFS_LMMSE(m,:) = berCOTFS_LMMSE(m,:)+errorStats_coded2.';  
        errorStats_coded1 = errorRate1(dataBits_in,dataBits_out1,1);          % 重置
        errorStats_coded2 = errorRate2(dataBits_in,dataBits_out2,1);
    end
    
end
%% 计算误码率
for i = 1:size(SNR,2)
    berCOTFS_MRC(i,1) = berCOTFS_MRC(i,2)/berCOTFS_MRC(i,3);
    berCOTFS_LMMSE(i,1) = berCOTFS_LMMSE(i,2)/berCOTFS_LMMSE(i,3);
end
%%     画图                  
plotGraphs1(berCOTFS_MRC, berCOTFS_LMMSE ,SNR_dB);
% save('mrc_lmmse.mat')
